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Abstract 

The MAP kinase cascade is a network of enzymatic reactions arranged 
in layers. In each layer occurs a multiple futile cycle of phosphorylations. 
The fully phosphorylated substrate then serves as an enzyme for the layer 
below. This papers focusses on the existence of parameters for which 
Hopf bifurcations occur and generate periodic orbits. Furthermore it is 
explained how geometric singular perturbation theory allows to generalize 
results from simple models to more complex ones. 


1 Introduction 

The MAP kinase cascade (mitogen-activated protein kinase cascade) is 
a pattern of chemical reactions encountered frequently in cell biology. 
The actual substances involved in the reactions vary from one example to 
another but certain features are always present. There are three proteins, 
generically denoted by MAPK, MAPKK and MAPKKK. For brevity we 
mostly use the notations K, KK and KKK for these. A phosphate group 
can be attached to KKK at a particular site and this causes it to become 
activated. It then catalyses the addition of phosphate groups to KK 
at two sites. An enzyme which catalyses a phosphorylation in this way 
is called a kinase and MAPKKK stands for MAPKK kinase. Similarly 
when KK has been phosphorylated at both sites it becomes activated and 
catalyses the addition of phosphate groups to K at two sites. There are 
also other enzymes (phosphatases) which remove phosphate groups from 
each of the three proteins. We refer to one of the proteins K, KK or KKK 
together with its phosphorylated forms and the reactions interconverting 
these forms as a layer of the cascade. We think of the layers corresponding 
to KKK, KK and K as being arranged from top to bottom and refer to 
them as the first, second and third layers. Since substances in one layer 
directly influence the reactions in the next layer information flows from 
the top to the bottom in this picture and that is the reason it is called a 
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cascade. The MAPK cascade can be represented in the following way: 


KKK 


KK 



( 1 ) 


KPP 


KP' ase 


KP' ase 


In this diagram, an enzymatic reaction X - 
ical reactions X-\-Et^X-E^Y-\-E, 


-s- Y represents the chem- 

where X ■ E \s a. complex of 


substrate X and enzyme E. In fact information can flow upwards through 
the cascade in a more indirect sense and this is described in more detail 
below. This backward flow has an essential influence on the dynamics of 
the system. 

If each of the reactions is modelled by a standard irreversible Michaelis- 
Menten scheme (see for instance [^) composed of three elementary re¬ 
actions with substrate, enzyme and substrate-enzyme complex and the 
elementary reactions are given mass action kinetics a system of ordinary 
differential equations arises as a model for the MAPK cascade. This 
will be referred to as the Huang-Ferrell model since it was introduced by 
those authors in m In that paper they studied mathematical properties 
of these equations and also compared solutions of the equations with the 
results of the experiments they had done on a MAPK cascade in extracts 
of oocytes (immature egg cells) of the frog Xenopus. The full system 
consists of 22 equations (for 8 substrates, 4 enzymes which do not occur 
as substrates and 10 complexes) and depends on 30 parameters (reaction 
constants). The procedure in [T7| was to fix all but one of the parameters 
and numerically determine stationary solutions of the system for different 
values of the remaining parameter. Then the value of the concentration of 
one of the substances at the stationary solution was studied as a function 
of the chosen parameter. This function exhibits the property of ultrasen¬ 
sitivity where the output is a sigmoid function of the input. 

The possibility of carrying out this procedure is dependent on the fact 
that there exist a stationary solution for each value of the parameters. 
Moreover the results will only be unambiguous if there is only one such 
solution for fixed values of the reaction constants and the total amounts of 
the three proteins and the four other enzymes. This issue is not addressed 
in m- That the answer to the uniqueness question is not obvious is made 
clear by the results of [22]. That paper was concerned with a system 
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which can be thought of as a single layer of the MAPK cascade with two 
phosphorylation steps. The ODE system in this case is known as the dual 
futile cycle m- The numerical and heuristic work in [22] indicated that 
this system exhibits bistability, i.e. that there are parameter values for 
which there exists more than one stable stationary solution. A rigorous 
proof that this is the case was given in m- 

A priori it is not ruled out that solutions of the Huang-Ferrell model 
might exhibit more complicated long-time behaviour than just conver¬ 
gence to a stationary solution. In fact, numerical and heuristic work in 
|26 | indicates that periodic solutions exist. There is also evidence suggest¬ 
ing that chaotic behaviour may occur 13S]- According to the investigations 
of [26] periodic solutions already occur in the system corresponding to the 
first two layers of the MAPK cascade. We refer to this system as the trun¬ 
cated Huang-Ferrell model. In what follows we will prove that for both the 
truncated and full Huang-Ferrell models there exist parameters for which 
there are periodic solutions. This means that the protein concentrations 
being modelled undergo sustained oscillations. 

Following the terminology of |14| we refer to the Huang-Ferrell model 
or an analogous system for a subset of the layers of the cascade as the 
MM-MA system (for Michaelis-Menten via mass action). Under suitable 
circumstances it is possible to derive a smaller system, the MM system 
(for Michaelis-Menten), via a quasistationary approximation. This is rel¬ 
atively simple to do for the dual futile cycle and the dynamics of the MM 
system in that case was studied in [24]. In that paper the authors found 
that bistability is already present in the MM system. The method used 
in m, which will be generalized here, is to first prove bistability for the 
MM system and then use the fact that the MM system is a limit of the 
MM-MA system in a suitable sense to obtain the corresponding result 
for the MM-MA system. For this we used that the stationary solutions 
are hyperbolic. The technique applied to treat the singular limit is ge¬ 
ometric singular perturbation theory (GSPT) [9]. In [16] an MM limit 
was defined for the truncated Huang-Ferrell model. The definition was 
inspired by ideas in [30] and [32]. In these papers it was pointed out that 
the phenomenenon of sequestration can lead to a flow of information from 
layers further down in the cascade to higher levels. For instance, if a lot 
of KKK*, the activated form of KKK, is bound to its substrates KK 
and the phosphorylated form KKP then not much of it is available as a 
substrate for the phosphatase which would otherwise convert it back to 
the inactivated form KKK. Related ideas are discussed in [6], where this 
effect is called retroactivity. It was mentioned in |16| that if it could be 
proved that this MM system admits a hyperbolic periodic solution then 
the truncated Huang-Ferrell system would also admit a periodic solution. 
It turned out to be difficult to obtain a proof of hyperbolicity and thus 
we will use a slightly different strategy in what follows. The basic idea is 
nevertheless still to first prove a result for the MM system and then lift 
some of the structure found to the MM-MA system. 

The paper is organized as follows. Section 2 introduces some basic 
notation and terminology and Section 3 proceeds to show that the MM 
system derived from the truncated Huang-Ferrell model admits a Hopf 
bifurcation for certain values of the parameters. It is then shown that 
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it can be concluded that the MM-MA system also admits a Hopf bifur¬ 
cation. It follows from the basic theorem of Hopf (cf. [TS]) that there 
are parameters for which the MM-MA system admits periodic solutions. 
These results are summed up in Theorem [T] Section 4 shows that the 
bistability in the dual futile cycle implies the presence of bistability in the 
truncated Huang-Ferrell model. In Section 5 the arguments of Section 3 
are extended so as to prove Theorem [2] which asserts the existence of peri¬ 
odic solutions for the full Huang-Ferrell model. In particular this involves 
the use of an MM system for the full MAPK cascade. Section 6 presents 
results on some variants of the models coming from the MAPK cascade. 
For a cascade consisting of two single phosphorylation loops it is proved 
that all solutions of the MM system converge to the same stationary solu¬ 
tion. In particular there are no sustained oscillations for that model. For 
a cascade consisting of a layer with two phosphorylations above a layer 
with one phosphorylation it is not clear whether sustained oscillations can 
occur. It is shown how this configuration can arise by modelling a system 
considered in [2^, thus motivating further study of this question. The 
final section indicates some directions in which the results of this paper 
might be extended. Some information on geometric singular perturbation 
theory is collected in an appendix. 


2 The basic equations 

In what follows we essentially adopt the notations of Huang and Ferrell 
m for the full cascade and then specialize them to the truncated system. 

The unknowns in the system are the concentrations 

[KKK], [KKK*], [El], [E 2 ], [KKK ■ Ei], [KKK* ■ E 2 ], 

[KK], [KKP], [KKPP], [KKP'ase], 

[KK ■ KKK*], [KKP ■ KKK*], [KKP ■ KKP'ase], [KKPP ■ KKP'ase[, 

[K], [KP], [KPP], [KP'ase[ 

[K ■ KKPP], [KP ■ KKPP], [KP ■ KP'ase], [KPP ■ KP'ase[. (2) 

The square brackets around a symbol indicate the concentration of the 
substance denoted by that symbol. Here Ei and E 2 are the kinase and 
the phosphatase in the first layer while KKP'ase and KP'ase are the 
phosphatases in the second and third layers. KKK* is the activated 
(i.e. phosphorylated) form of KKK. A P occurring in the name of a 
protein indicates a phosphate group. For instance KKPP is the doubly 
phosphorylated form of KK. Finally, the substrate-enzyme complexes 
are denoted by the symbols for the two substances separated by a centred 
dot. The full set of 22 evolution equations for these quantities will not be 
reproduced here. They can be found in m- More precisely, 18 of them 
can be found there and the other four, those for the concentrations of the 
free enzymes, are easily derived from those. There are seven conserved 
quantities, three corresponding to the total amounts of the three substrate 
proteins in all phosphorylation states and four corresponding to the total 
amounts of the four enzymes. These conservation laws can be substituted 
into the evolution equations so as to reduce their number. In m the four 
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conserved quantities associated to the enzymes are used in this way while 
the others are not. Thus 18 evolution equations remain. The reaction 
constants for the elementary reactions in the Michaelis-Menten description 
are denoted by hi for the formation of the complex, ki for the liberation of 
product and di for the release of the substrate from the complex. These 
notations are as in m except that at and di have been replaced by hi and 
di since ai and di will be used for a different purpose later. All reaction 
constants are assumed positive and since the unknowns in the equations 
are concentrations the solutions of interest are those which are positive, 
i.e. all their components are positive. 

The truncated Huang-Ferrell model is obtained by setting the concen¬ 
trations of all the unknowns in the full model involving the substance K 
and that of the enzyme KP'ase (i.e. the quantities in the last two lines 
of Q) to zero and discarding the evolution equations for those quanti¬ 
ties. For the rest of this section we concentrate on the truncated model. 

This serves, in particular, to introduce some of the main techniques of the 
paper in a context simpler than that of the full cascade. The results for 
the truncated cascade also constitute a central component of the proof 
of the theorems for the full cascade. As explained in [16] the passage to 
the limiting MM system can be carried out by scaling the variables in a 
certain way with powers of a parameter e and then letting e tend to zero. 

First it is important to introduce a new variable to replace [KKK*\. This 
is 

KKK = [KKK*] + [KK ■ KKK*] -h [KKP ■ KKK*]. (3) 

When the Huang-Ferrell model is written in terms of KKK the equa¬ 
tions of the truncated model are a subset of the equations for the full 
model. The scaling is done as follows. The quantities in the evolution 
equations are replaced by new quantities defined in the following way. 

The concentrations of compounds not containing KKK or the enzymes 
Ei or KKP'ase are not rescaled. These are the first three quantities in 
the second line of For the concentrations of compounds which contain 
KKK or KKP'ase but not the Ei the new quantity is e~^ times the old 
one. These are the first two quantities in the first line, the last quantity 
in the second line and all the quantities in the third line of m- For the 
quantities involving Ei and E^ the new quantity is times the old one. 

These are the last four quantities in the first line of The reaction con¬ 

stants hi and 02 are multiplied by e to get new quantities. All the other 
reaction constants are left unchanged. In addition a new time coordinate 
is introduced as e times the old one. To avoid complicating the notation 
the rescaled quantities will be denoted in the same way as the original 
ones. There results the following system, which was given in a different 
notation in |16| . 

^^(KKK) = -h2[KKK*][E2\ + MKKK* ■ E 2 ] + ki[KKK ■ Ei], (4) 

^[KK] = -h 3 [KK][KKK*] + dslKK ■ KKK*] + ki[KKP ■ KKP'ase], (5) 

^[KKPP] = -h6[KKPP][KKP'ase] + deiKKPP ■ KKP'ase] + k^lKKP ■ KKK((I} 
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^ di[KKK][Ei] - (di + ki)[KKK ■ Ei], (7) 

^ d 2 [KKK*][E 2 ] - {d 2 + k 2 )[KKK* ■ E 2 ], (8) 

dt^^ ^ "" ^3[KK][KKK*] - {ds + k3)[KK ■ KKK% (9) 

KKK*] ^ a3[KKP][KKK*] - (dg + k3)[KKP ■ KKK*], (10) 


^ d[KKP ■ KKP ase] ^ a4KKP][KKP'ase] - {d^ + k4)\KKP ■ KKP'ase], (11) 
^ d\KKPP -KKP ase\ ^ KPP][KKP'ase] - {d& + ke,)[KKPP ■ KKP'ase]. (12) 

These equations correspond to a subset of the equations given in [T^ ■ To 
get a closed system we must express some of the quantities on the right 
hand side of the equations in terms of those on the left hand side using the 
definition of KKK and the conservation laws. The necessary equations 


are 

[KKK*] = KKK - ]KK ■ KKK*] - [KKP ■ KKK*], (13) 

[KKK] = KKKtot - liKK + 0(e), (14) 

[KKP] = KKtot - [KK] - [KKPP] + 0{e), (15) 

[Ei] = Ei,tot-[KKK-Ei], (16) 

[E 2 ] = E 2 ,tot - [KKK* ■ E 2 ], (17) 


[KKP'ase] = {KKP'a.se)tot - ]KKP ■ KKP'ase] - ]KKPP ■ KKP'c{A^) 


This system of equations extends smoothly to e = 0. The terms written 
as 0{e) are sums of concentrations of complexes. They are not written 
explicitly since they vanish for e = 0 and thus make no contribution to the 
limiting equations. When e = 0 the equations ([7l)-(fT^ become algebraic 
equations. Substituting these into the other evolution equations gives 
a set of three evolution equations, the MM system, which will now be 
computed. First the concentrations of the complexes can be expressed in 
terms of those of the substrates and the free enzymes. 


[KKK • Si] = 


[KKK* ■ E 2 ] = 
[KK ■ KKK*] = 


di[KKK][Ei] 
di + ki 

a 2 [KKK*][E 2 ] 


d2 + k2 

a 3 [KK][KKK* 


[KKP ■ KKK*] = 


ds + /l3 

a4KKP][KKK*] 


[KKP ■ KKP'ase] = 


ds + fcs 

a4KKP][KKP'ase] 


[KKPP ■ KKP'ase] = 


d4 “t” 

ae[KKPP][KKP'ase] 


dg “t“ /bg 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 


It is convenient to introduce the Michaelis constants Km,i = . Now 

the concentrations of the free enzymes can be expressed in terms of the 
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total concentrations of the enzymes. Note that 


+ K-\[KKK]), (25) 

£;2.tot = [£;2](1 + K-]^[KKK*]), (26) 

WKK = [KKK*]{1 + K-^s[KK] + K-^^[KKP]), (27) 

{KKP'ase)tot = [KKP'ase]{l + K-\[KKP] + Jf-ig [Tfif PP])28) 


Substituting this into the equations 


^(KKK) = -^[KKK^][E2] + -^[KKK][E^l (29) 

at Am,2 Am,l 

^ [A'P:] [PTP'P'* ] + [Tt'PTP] [P'P'P'ose], (30) 

dt Am,3 Am,4 

^[P'P'PP] = --^[7s'P'PP][p:P'P'ase] + -;:^[KKP][KKi&l]) 
dt Km,6 Km,5 


gives 


d ^ k.K-\E.MKKK] 

dt'' 1 + K-\[KKK] 

A )2 Ayj^ 2 A2, t ot KKK 

~ 1 + K-:,KKK + K-‘:,[KK\ + K-^,[KKP] ’ 

I = k3K-],KKK[KK] 

l + i^-;3[7^i^]+i^-;5[i^7^P] 


l + K-^,[KKP] + K-^,[KKPPy 
d , , k5K;;yj<KK[KKP] 

dt' l + K;y^^[KK] + K;;yyK 


iKKP'ase)tot [KKPP] 

l + K^y[KKP] + K^\[KKPP] 


(32) 


(33) 


(34) 


Equations (EU-dMI are the Michaelis-Menten system. For convenience we 
introduce the notations ai = fci+27f^\+2 for i = 1, 3, fli = ki+ 2 K^\^ 2 {d^diP'ase)tot 
for i = 2, 4, Ci = kiK^dEi^tot for i = 1, 2 and di = for i = 1, 2. Fur¬ 
thermore we replace [A'A'P], [KKK] by their expressions in the chosen 
phase variables KKK, \KK], [A'A'PPJand the parameters of the system. 

Assume that the coefficients K^^ are equal to a common quantity b\ 
for 3 < i < 6. For any choice of the parameters ai, bi, a and di there 
exists a choice of the parameters and total quantities of the enzymes for 
the MM-MA system which gives rise to them. This choice can be made 
so as to depend smoothly on {ai,b\,Ci,di). 

In this notation, the Michaelis-Menten system for the truncated cas¬ 
cade is the following: 


-(Tan<) = cijKKKtot-KKK) 
dP 1 + di {KKKtot - ICKK) 
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C2KKK 

1 + d2KKK + bi {KKtot - [KKPP] ’ 


( 35 ) 


d a{KKK[KK] 

de ^ l + hi{KKt^t-\KKPP]) 
a2{KKtot - [KK] - [KKPP]) 
l + bi{KKtot-[KK]) 

^IKKPP] = ^3Tnn<{KKtot - [KK] - [KKPP]) 
dP ^ l + hpKKt^t-[KKPP] 

a^KKPP] 

~ I+ bi{KKtot - [KKPP])' 


3 Analysis of the truncated system 


The aim of this section is to prove that there are values of the parameters 
and the total amounts of substrates and enzymes in the truncated Huang- 
Ferrell system for which there exist periodic solutions. The first step is to 
determine explicit stationary solutions of the system (I551)-(I571) for certain 
values of the total amounts. The existence of periodic solutions of the MM 
system is obtained by showing that a Hopf bifurcation occurs at these sta¬ 
tionary solutions. It is then shown that there is also a Hopf bifurcation, 
and hence periodic solutions, for the full truncated Huang-Ferrell system. 
Suppose that all parameters Ui, 6 i, C 2 and d; have been fixed. Consider 
stationary solutions of (I35I) - (I37I) which satisfy [KK\ = [KKPP\. Then 
a{KKK\KK] = a 2 [KKP] and a 3 KKK[KKP] = aA\KKPP] as a conse¬ 
quence of ([551) and (1371) . In particular the product of these two equations 
_ 2 _ 

gives KKK = and so KKK is determined. Next it is possible 

to determine and \KKP] as functions of KKtat by means of the 

formulae 


[KK] 


KKtot 
2 + qi’ 


[KKP] 


q±K Ktot 

2 qi 


(38) 


where qi — y . The first equation (1351) has not been used yet and 
produces a constraint on the parameters. This is the price to pay for 
simplifying some computations by having an equilibrium whose two last 
coordinates are identical. Set Q — , , ^ It follows 

l-\-d2KKK-\-bi {KKtot — [i^K\) 

from the equation dKKK/dt — 0 that , = Q- This is 

' l+di{KKKtot—KKK) ^ 

equivalent to KKKtot = KKK-\ -^- 7=7 provided the denominator in the 

Cl fii V 

last term is positive. Choose ci so that <; \ Then the denominator 
is positive and this expression can be used to define KKKtot- Once this 
has been done all the conditions for stationary solutions are satisfied. Let 



























US summarize this information: 


coordinates of the equilibrium 


constraints on the parameters 


KKK = 

' ' aia 3 


[KK] = \KKPP] = 

£2dl < I 


KKKtot = 

where Q = 


<t2Q4 _|_ w 
aias Cl—diQ ’ 


c^KKK 


l+d2KKK+bi{KKtot-[KK]) 

(39) 


Next the linearization of ISSli-dszl) will be considered. It turns out 
that the signs of all entries in the derivative of the right hand side of the 
system are independent of the values of the concentrations and can be 
determined. One is zero, one is positive and all the rest are negative. 
The characteristic polynomial of the linearization is the determinant of a 
matrix of the form 


A + a 0 c 
d X + e f 
—g h A + i 


(40) 


where all parameters are positive. Note that —e, —/, —h and —i only 
differ from the elements of the matrix (46) of [TS] by the multiplicative 
factor N~'^ and the fact that ai and as are replaced by a\KKK and 
asKKK, respectively. As in that case the parameters aiKKK and 04 
can be eliminated. Let M 2 be the submatrix of the matrix (|iD|l with A = 0 
consisting of the second and third rows and columns. The determinant of 
the linearization is given by 


— {a\M 2 \ + c{dh + eg)}. 


(41) 


Since we are seeking parameters for which a Hopf bifurcation producing 
stable periodic orbits takes place, the three eigenvalues of the linearization 
at the equilibrium (I39II should be two imaginary eigenvalues Piuio, cuq yf 0, 
and a third real negative eigenvalue. In particular the determinant of the 
linearization and its trace should be both strictly negative. The sign 
of I M 21 is unclear, but the linearization has a negative determinant if 
IM 2 I is small enough. Define u = hiKKtot and v = bi[KK]. Then the 
determinant of M 2 is given, up to a positive multiplicative constant, by 

— (u — 3u + l)(u — 2v)u{u — u + 1). (42) 


In fact all factors except the first are positive for all allowed values of the 
variables. Thus the determinant is negative iff u — 3w + 1 >0. This is 
equivalent to the condition ei — fh < 0. In terms of the parameters of 
the reduced system u — 3u + 1 = ^ biKKtot + 1)^. It follows that 

provided gi < 1 the value of KKtot can be varied so that the determinant 
of M 2 passes from being positive to being negative. 

Now consider the determinant of the matrix in the linearized system. 
It is equal to —a(ei — fh) — c{dh + eg). The second term is automatically 
negative and so if ei — //i > 0 the determinant of the full matrix is neg¬ 
ative. On the other hand if ei — fh < 0 the sign of the determinant is 
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not clear. The trace is negative for all parameter values. In a deforma¬ 
tion as above the determinant of the three-dimensional matrix is negative 
when the determinant of M 2 becomes zero. It stays negative at least in a 
small neighbourhood of that point. Thus there are parameters for which 
both the determinant of M 2 and the determinant of the three-dimensional 
matrix are negative. The determinant of PI|) is 

+ {a + e + i)\^ + {ae + ai + cg + ei —fh)X + {aei — afh + cdh + cge). (43) 

Call it P 3 (A). The eigenvalue conditions for a Hopf bifurcation are that 
two eigenvalues are purely imaginary and non-zero and that the third 
is non-zero. If this is true the trace is given by the real eigenvalue and 
thus determines its sign. In the example here we see that at a Hopf 
bifurcation the real eigenvalue must be negative. The determinant is then 
also negative. Introduce the notation P 3 (A) = A® -|- A 2 X^ -|- AiA -I- Ao- In a 
deformation as above all Ai are positive near the point where ei — fh — 0. 
The Routh-Hurwitz theorem m implies that all roots of this polynomial 
have negative real parts if and only if Ai > 0 for all i and H 3 = A 2 A 1 — 
Ao > 0. The Routh-Hurwitz coefficient Hs is given by 

H 3 = (a -\- e -\- i){ae + ai + eg + ei — fh) — aei + afh — cdh — ege (44) 
which simplifies to 

H 3 = a{ae + ai + eg) -|- (e -|- i){ae -I- a* -|- ei) + cgi — [cd -I- (e -I- i)f]h. (45) 

Let S 3 be the region in the space of three by three matrices defined by 
the condition that all eigenvalues have negative real parts. At any point 
on the boundary of S 3 where the determinant is non-zero the eigenvalue 
configuration characteristic of the Hopf bifurcation must occur. In other 
words, two eigenvalues are imaginary and non-zero and the other is neg¬ 
ative. 

In the Michaelis-Menten system for the truncated MARK cascade we 
have already exhibited a class of stationary solutions and shown that 
there is a choice of the parameters for which the matrix defined by the 
linearization of this system has negative determinant and the determinant 
of M 2 is negative. Call this point in parameter space Z. We recall that 
at Z, the constraints that have to be fulfilled are the following: 

Cid2 

KKKtot 

gi 

KKtot 

where <5 > 0 is small. At Z the linearization has negative trace. The 
coefficients will now be rescaled in a particular way. Let di be the values 
of the parameters ai at the starting point and define ai{L) = diL. The 
other parameters are kept fixed. Under this rescaling a stationary solution 
remains a stationary solution. The determinant of the linearization is 
rescaled by and so, in particular, has a hxed sign. The constants d- 
i in the matrix of the linearization are scaled by a factor L while the 


< 


a2Q4 

<11113 


£ 1£4 ^ I 
< 12<13 

3iA2.^+S, 


(46) 
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coefficients a and c are left unchanged. After rescaling the condition that 
the quantity H 3 is negative becomes 

L^( fh — ei){e + i) > La{ae + ai + eg) + L^(e + i){ae + ai) + L^cgi — L^cdh. 

(47) 

If we start the rescaling from Z then for L sufficiently large this condition 
will be satisfied. On the other hand, for L sufficiently small and positive 
the opposite inequality will be satisfied. In other words, in this family the 
quantity Hz starts positive for L small and positive and becomes negative 
at some large value of L. It must pass from being positive to being zero at 
some point Lq. In fact it must have a minimum for some negative value 
of L and a maximum for some positive value of L since it vanishes at 
the origin and its derivative there is positive. In fact the value Lq is the 
unique positive root of a quadratic polynomial. At the unique positive 
value I/O of L for which it is zero this function has a negative derivative. 
It follows from the facts that A 2 > 0 and Aq > 0 that Ai cannot approach 
zero before Hz reaches zero. Thus the family cannot leave Sz before Hz 
reaches zero. This family passes through a point where the eigenvalues 
satisfy the conditions for a Hopf bifurcation. When a three-dimensional 
matrix has one negative real eigenvalue and another eigenvalue ^2 
which is not real the quantity Hz is equal to 

— Rep2[(ReAi2)^ -b 2(Imp2)^ -b + Re^2)^]. (48) 

From this we can see that if a one-parameter family of matrices is such that 
the quantity Hz passes through zero at some parameter value with non¬ 
zero velocity then the real part of the complex eigenvalue passes through 
zero at that point with non-zero derivative. 

The situation which has just been described is as follows. We have 
a family of coefficients in the Michaelis-Menten system depending on the 
parameter L. There is a corresponding family of stationary solutions. We 
consider the linearization of the system at that point. We can suppose 
that for all values of L considered it has one negative real eigenvalue pi 
and one eigenvalue p 2 which is never real. Moreover at L = Lq the eigen¬ 
value p 2 is imaginary and its real part moves through the imaginary axis 
with non-zero velocity as L is varied through Lq. In other words, fJ. 2 {L) 
defines a curve in the complex plane which passes through the imaginary 
axis transversely for L = Lq. In particular the linearization is always 
invertible. It was shown in m that the Michaelis-Menten system ( 1551 ) - 
(1371) can be embedded into the MM-MA system ([i |) - (ll2l) in such a way as 
to allow GSPT to be applied. In the terminology of [16] the transverse 
eigenvalues all have negative real parts. More information on these mat¬ 
ters can be found in the appendix. Above we introduced a parameter L 
into the MM system by letting the coefficients depend on L in a certain 
way. This can be extended to the MM-MA system by making a smooth 
choice of corresponding coefficients for that system. After doing this the 
system on the slow manifold depends on the two parameters e and L. At 
the point where the Hopf bifurcation occurs for e = 0 and L = Lq the 
derivative of the right hand side of the system with respect to the un¬ 
knowns is invertible. Hence, by the implicit function theorem there is a 
unique nearby stationary solution depending smoothly on L and e. Thus 
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the family of stationary solutions for e = 0 can be continued smoothly to a 
curve of stationary solutions in the slow manifold for small positive values 
of e. For each e there are corresponding eigenvalues 
and fl 2 {L,e) = of the linearization at the stationary solution. For 

e small the curve ^ 2 ,e intersects the imaginary axis transversely. Now the 
transverse eigenvalues are negative and it follows that the family of solu¬ 
tions of the MM-MA system undergoes a Hopf bifurcation. Using centre 
manifold theory the dynamics can be reduced to a two-dimensional situa¬ 
tion for both the MM-MA system for the truncated Huang-Ferrell model 
and the corresponding MM system. Applying the theorem of Hopf in the 
form given in m, Chapter VIII, Theorem 1.3 shows that in both cases 
there exists a one-parameter family of periodic solutions. The following 
theorem has been proved. 

Theorem 1 There exist positive parameter values at, bi, a and di such 
that the MM system for the truncated Huang-Ferrell model with these pa¬ 
rameter values has a positive periodic solution. There exist positive param¬ 
eter values di, ki and di and positive values of the total amounts 
E 2 ,tot, {KKP'ase)tot, KKKtot and KKtot such that the MM-MA system 
for the truncated Huang-Ferrell model with these parameter values and 
these values of the total amounts has a positive periodic solution. 


4 Bistability 

In [16] it was proved that there are parameter values for which the dual 
futile cycle has more than one stable stationary solution. It will now be 
shown that this property is inherited by the truncated MAPK cascade. 
For this the coefficients in the MM system for the truncated MAPK cas¬ 
cade will be rescaled. Define ci = C 2 = e~^C 2 and d 2 = e~^d 2 . 

Then the equations for the quantities with hat are in the standard form 
for singular perturbation theory. For fixed values of KKK the equations 
for \KK] and [KKPP] are those for the dual futile cycle and we know 
that there are parameters for which there exist two stable hyperbolic sta¬ 
tionary solutions. The evolution equation for KKK is 

^ ci(AAAtot - 

l + di{KKKtot-lMn<) 
_ 02 AAA _ 

e -b d2KKK -b e6([AA] + [AAP])' 

For e = 0 this reduces to 

di (^ -1 j {KKKtot - Taac) = 1. 


(49) 


(50) 


Interestingly this last equation does not depend on [AA] or [AAP]. Pro¬ 
vided the quotient of reaction constants occurring in this equation, which 
is equal to is greater than one then a unique stationary solution is 

determined by solving this equation for KKKtot in terms of KKK. It is 
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an asymptotically stable and hyperbolic solntion of the eqnation 


( 51 ) 

l + di{KKKtot-KKK) d2 

Thus the transverse eigenvalue in the sense of GSPT is negative and that 
theory can be applied. It follows that the MM system for the truncated 
MAPK cascade exhibits bistability. This implies a corresponding state¬ 
ment for the MM-MA system. It can be shown in a similar way that both 
the MM and MM-MA systems contain saddle points for these values of 
the parameters. 

It can also be shown that there are parameters for which there is a sta¬ 
tionary value where the eigenvalues are real with the sign pattern ( —, —, 0). 
There ei — fh < 0. Scaling with the parameter L starting at this point 
gives a family of coefficients for which the determinant of the linearization 
is always zero. Along this family the quantity H 3 goes from positive to 
negative. When it crosses zero A 2 = 0 and there is a second zero eigen¬ 
value. If the kernel of the matrix at that point were two-dimensional it 
would have to intersect the subspace spanned by the first two coordinates 
in a one-dimensional subspace. It would follow that this subspace was 
spanned by the vector with components (0, 1,0). This gives a contradic¬ 
tion. Thus the Jordan form of the matrix must be non-diagonalizable 
and the algebraic conditions on the linearization for a Bogdanov-Takens 
bifurcation are satisfied at that point. In other words, the condition BT.O 
of [ 1 ^ holds. Since, on the other hand, no information has been obtained 
on the genericity conditions BT.1-BT.3 this does not by itself give useful 
information on the dynamical properties of the solutions. 


5 The full Huang-Ferrell system 

The aim is now to reduce the full Huang-Ferrell system in a way similar to 
that done above for the truncated system. The first step is to introduce 
the quantity 

liK = [KKPP] + [K ■ KKPP] + [KP ■ KKPP]. (52) 

Consider the system of evolution equations for KKK, [KK], KK, [K\, 
[KPP] and the substrate-enzyme complexes and rescale the unknowns. In 
this case the quantities not containing the Ei, KKK, KK, KKP'ase or 
KP'ase are not rescaled. These are the first three quantities in the fourth 
line of ©■ Quantities which contain KK or KP'ase but not Ei, KKK 
or KKP'ase are rescaled by e~^. These are the first three quantities in 
the second line, the last quantity in the fourth line and all quantities in 
the fifth line of ((2]). Quantities which contain KKK or KKP'ase but not 
Ei are rescaled by e~^. These are the first two quantities in the first line, 
the last quantity in the second line and all quantities in the third line of 
©. Quantities which contain the Ei are rescaled by e These are the 
last four quantities in the first line of @. The reaction constants hi and 
0,2 are multiplied by to get new quantities while 0,3, 0,4, dz and he are 
multiplied by e. A new time coordinate is introduced as e times the old 
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one. These scalings have been chosen so that the new eqnations for the free 
substrates are independent of e and the new equations for the substrate- 
enzyme complexes have a factor e in front of the time derivatives, just as 
in the case of the truncated system. The resulting equations extend ©- 
dH. The first two rescaled equations are unchanged and the third only 
differs from the corresponding equation for the truncated system in that 
[KKPP] is replaced by KK on the left hand side. The list of expressions 
for the concentrations of the complexes can be extended as follows 

[K.KKPP] = ^-im^, 

lKP.KP'ase]=^-^hM^, 
ds ^8 

[KP.KKPP] = ^-^im^^^, 

dg “h kg 

[KPP ■ KP’ase] = PP][KP'ase] 

dio -|- fcio 

There are the following additional equations for the total amounts of en¬ 
zymes 

-KK = [KKPP](1 + (57) 

{KP'ase)tot = [KP’ase]{l + K;;^\[KP] + (58) 

The equations (p9l) -([3T|) can be taken over to the full model except that in 
m the quantity [KKPP] should be replaced by KK on the left hand side 
and that in order to obtain a closed system [KKPP] must be substituted 
for in terms of KK on the right hand side. The following equations also 
hold: 

^ [7^1 = - [K] [KKPP] + ]KP] [KP'ase ], (59) 

at i\m,7 Am,8 

^[P'PP] = --^^]KPP]]KP'ase] + -^[KP]]KKPP].{m) 

dt Km,10 Km,9 


(53) 

(54) 

(55) 

(56) 


The evolution equations for the MM system are 

d ^ k.K-J,E.MKKK] 

dt'' 1 + K-\]KKK] 

k2K-^2E2,tojn{K 

1 + K-^,-KKK + K-^,]KK] + K-],]KKP] ’ 
^ = k3K-%-KKK[KK] 

^ + k-%[kk]+k-^,[kkp] 

k4K-^^{KKP'ase)tot]KKP] 

+ 1 + K-^^]KKP] + K-],]KKPP] ’ 

^ k3K-\-KKK]KKP] 

1 + K-^^]KK] + K-^,]KKP] 


(61) 


(62) 
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keK-]e{KKP'ase)tot[KKPP] 

~ 1 + K-^4[KKP] + K-^,[KKPP] ’ 

^ k7K-],in<[K] 

di 1 + K-^AK]+K-^9[KP] 


+ 


k8K-^siKP'ase)tot[KP] 

^ + k;;.^[kp] + k-\,ikpp]’ 


kioK-\o{KP'ase)tot[KPP] 

^ + K;;.]slKP] + K-\o[KPPy 


( 63 ) 


(64) 


(65) 


Assume that the coefhcients K^\ are equal to a common quantity 62 for 
7 < i < 10. Extend the definition of at by defining it to be fci+ 2 A '“\_|_2 
for i = 5, 7 and fci+ 2 A“f_,_ 2 (AP'ase)tot for i = 6 , 8 . As in the case of 
the truncated system there exists a smooth choice of the parameters and 
total quantities of the enzymes for the MM-MA system which give rise to 
any choice of the parameters at, bi, a and di. 

We now have a system describing the full cascade which is in the stan¬ 
dard form of GSPT. To profit from this it is necessary to examine the 
transverse eigenvalues. These are the eigenvalues of the matrix which 
is the derivative of right hand side of the evolution equations for the 
enzyme-substrate complexes for fixed values of the substrate concentra¬ 
tions. There are ten complexes and this is a ten by ten matrix with 
components Lij. Let the complexes be numbered in the order they are 
listed in ©. Each complex which does not share an enzyme contributes 
a diagonal element to the matrix. These are the components Ln and 
L 22 and are negative. Each pair of complexes which share an enzyme 
contributes a two by two submatrix on the diagonal. They are the Lij 
with 2k — 1 < i, j < 2k for 2 < k < 5. The eigenvalues of each of these 
submatrices on the diagonal have negative real parts. The calculation is 
essentially the same as that done for the case of the truncated system in 
m- It remains to examine the effect of the non-zero elements of the ma¬ 
trix which do not belong to any of these blocks. These are L 23 , I/24, Lez 
and 1/68- The elements Ln and L 22 are alone in their columns. Thus the 
calculation of the eigenvalues reduces to that of the submatrix obtained 
by discarding the first and second rows and columns. Then the subma¬ 
trices for k = 1 and fc = 4 occur as direct sums with other matrices. 
Thus determining the eigenvalues reduces to doing so for the submatrix 
defined hy 5 < i, j < 8. This submatrix is block upper triangular and so 
its eigenvalues are the eigenvalues of the submatrices for k — 2 and fc = 3. 
Combining these facts shows that all transverse eigenvalues have negative 
real parts. 

Consider a stationary solution of the system dSIli-dell) which satisfies 
[KK] — [KKPP] and [K\ = [KPP\. Explicit stationary solutions can 
be found in a way similar to what was done for the truncated system. 
It follows from equations (IMl) and (1651) that a 3 KK[K] = a8[KP\ and 
a 7 'KK[KP] = aslKPP], Hence 'KK'^ = The quantities [K\ and 
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[KP] are determined by 


[K] 


^tot riyp] _ Q2P-tot 

2 + 92’ ^ ‘ ~ 2 + q 2 


( 66 ) 


where q 2 = y The expression obtained for KKK in the case of the 
truncated system remains valid for the full system while in those for [KK] 
and [KKP] are modified to 

[KK] = [KKP] = (67) 

2 + qi 2 + qi 

where qi — qi + b 2 {Ktot — [7f])- In addition we have the relation 


KKtot 


{2 + qi)KK 
1 + &2(7ftot - [7f]) ’ 


( 68 ) 


so that KKtot is determined. Stationary solutions for the MM system 
for the full cascade can be determined as follows. Fix the parameters 
fli, hi, C 2 , di. Then if ci is chosen sufficiently large the concentrations 
for the stationary solution can be expressed in terms of Ktot and KKtot- 
Note that if the parameters tti, 1 < i < 8 are varied in such a way that 
rk = ci 2 k+ilO‘ 2 k remains unchanged for 1 < fc < 4 then qi and q 2 do not 
change and the stationary solution is preserved. 

Expressing the evolution equations in terms of the quantities at, hi, a 
and di and using the conservation laws gives 


_d ^ ct{KKKtot-innc) 

dK l + dt [KKKtot - 'KKK) 

C2KKK 

~ 1 + d2'KKK + &i [KKtot - Tn<) ’ 

— \KK] = ajniK[KK] 

dV l + bi{KKtot-in{) 

^ _ a2{KKtot - [KK]-in<) _ 

1 + 61 [KKtot - [KK] - ’ 

d — aaTOO?{KKtot - [KK] - ICK) 

dK ' l + bt{KKtot-[KK]) 

_ aJ<K _ 

b{KK + (1 + bt{KKtot - [KK] - 'KK )){1 + h2{Ktot 
d aii'KK[K] 

dP ^ l + h2{Ktot-[KPP]) 
ae{Ktot - [K] - [KPP]) 

I + b2{Ktot - [K]) ’ 

Ir^ppi ^ aj'KKiKtot - [K] - [KPP]) 
dP ^ ^ l + h2{Ktot-[KPP]) 

as[KPP] 

1 + h2{Ktot - [K]) • 


(69) 


(70) 



(72) 


(73) 
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A relation will now be established between this system and the system for 
the truncated model by doing some rescaling. Replace [K], [KPP] and 
Ktot by t[K], t[KPP] and eAtot, respectively. Replace at by t~^ai for 
5 < i < 8. In the limit the first three equations are independent of \K\ and 
[KPP] and are just the equations of the truncated system analysed in the 
last section. The limit is in the form appropriate for applying GSPT. The 
linearization of the system for [K] and [KPP] is independent of those two 
variables. Its trace and determinant are —{as + ar)KK — (ae + as) < 0 and 
a^a-jKK +(a 5 + a6)a8 > 0, respectively and so the transverse eigenvalues 
have negative real parts. A parameter L can be introduced in the system 
for the full cascade in the same way as was done for the truncated cascade. 
It follows that the presence of a Hopf bifurcation in the MM system for 
the truncated model implies that of Hopf bifurcation in the MM system 
for the full cascade. A parameter L can also be introduced in the MM-MA 
system for the full cascade, implying the existence of a Hopf bifurcation 
for that system, i.e. the original system of Huang and Ferrell. Thus the 
following theorem has been proved. 

Theorem 2 There exist positive parameter values at, bi, a and di such 
that the MM system for the full Huang-Ferrell model with these parame¬ 
ter values has a positive periodic solution. There exist positive parameter 
values hi, ki and di and positive values of the total amounts i? 2 ,tot, 

{KKP'ase)tot, {KKP'ase)tot, KKKtot, KKtot and Ktot such that the 
MM-MA system for the full Huang-Ferrell model with these parameter 
values and these values of the total amounts has a positive periodic solu¬ 
tion. 


6 Further examples 

This section is concerned with some examples which are variations on 
those coming from the MARK cascade. The first is a system which is 
similar to the truncated Huang-Ferrell system except that the second layer 
of the cascade only has one phosphorylation site. In other words, we 
discuss now the following cascade: 


El 


KKK KKK* 



KKK" 


(74) 


KK KKP 

KKP'ase 

This minimal cascade and its generalization to several layers of simple 
phosphorylation loops have been considered in [30] where it is remarked 
that damped oscillations may occur in a system of this type. It will now 
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be shown how the features observed for the MAPK cascade change in the 
case of a cascade of two simple phosphorylation loops. The variables are 


[KKK], [KKK*], [El], [E2], [KKK ■ Ei], [KKK* ■ E2], 

[KK], [KKP], [KKP'ase], [KK ■ KKK*], ]KKP ■ KKP'asep^) 

We introduce the variable KKK = ]KKK*] + [KK ■ KKK*] in analogy 
to what was done in the previous examples. The variables are rescaled as 
in the truncated Huang-Ferrell model. This means that [KK] and [KKP] 
are not rescaled, for [KKK], [KKK*], [KKP'ase], [KK ■ KKK*] and 
[KKP ■ KKP'ase] the new quantity is e~^ times the old one while for 
[^ 2 ], [KKK ■ E\[ and [KKK* ■ £ 2 ] the new quantity is times 
the old one. The reaction constants fii and 0,2 are multiplied by e to get 
new quantities. A new time coordinate is introduced as e times the old 
one. There results the following system 

^^(KKK) = -a2[KKK*][E2] + d2[KKK* ■ E2] + ki[KKK ■ Ei], (76) 

^[KK] = -a3[KK][KKK*[ + d3[KK ■ KKK*] + ki[KKP ■ KKP'ase]p 7 ) 

= di[KKK][Ei] - (di + ki)[KKK ■ Ei], (78) 

^dp<J<K_^ = d 2 [KKK*\[E 2 \ - {d 2 + k 2 )[KKK* ■ £ 2 ], (79) 

^d[KI<_KK^ = a3[KK][KKK*] - {ds + k3)[KK ■ KKK*], (80) 

^ d[KKP ■ KKP ase] ^ ai[KKP][KKP'ase] - {di + k4)[KKP ■ AAP'^M]) 

To get a closed system the following relations must be used. 


[KKK*] = KKK - [KK ■ KKK*], (82) 

[KKK] = KKKtot - TCKK + 0(e), (83) 

[KKP] = KKtot - [KK] + 0{e), (84) 

lEi]^Ei,tot-[KKK-El], (85) 

[E 2 ] = P 2 .tot - [KKK* ■ P 2 I, (86) 

[KKP'ase] = {KKP'ase)t„t - [KKP ■ KKP'ase]. (87) 


For e = 0 the following MM system is obtained 

^ ^ fciA-^iPi,tot[AAA] 

dP 1 + K-\[KKK] 

fc2A-;,P2 ,tot KKK 

~i + lp^^WWWTlp^^{KKy 
- k3K-]pnni[KK] 
dt^ ^ f + A-^AA] 

^ fc4A-i4(AAP'ase)tot[AAP] 

+ l + A-i^lAAP] • 


( 88 ) 


( 89 ) 
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The system (I88I) - (I89I) has the property that the trace of the derivative of 
the right hand side is always negative. Thus by the Dulac criterion this 
system admits no periodic solutions. It was shown in [S] that in this case 
the MM-MA system has a unique stationary solution for given values of 
the parameters and this means that the same is true for the MM system. 
On the boundary of the region of positive concentrations the vector field 
points inwards and all solutions are bounded due to the conservation laws. 
Putting all these facts together, it follows from Poincare-Bendixson theory 
that the stationary solution of the MM system is globally asymptotically 
stable. It is not clear that stability might not be lost for general values of 
the parameters in the MM-MA system. 

The next example is an in vitro model, introduced by Prabakaran, 
Gunawardena and Sontag [^, of the MAPK cascade consisting of the 
proteins Raf, MEK and ERK. The model system is simplified compared 
to the original biological system in two ways. The protein Raf is consti- 
tutively active. This corresponds to taking a fixed value of KKK in the 
model. ERK is mutated so that it can only be phosphorylated once, on 
tyrosine and not on threonine. The role of the phosphatase KKP'ase is 
played by PP2A (protein phosphatase 2A) and that of KP'ase by PTP 
(protein tyrosine phosphatase). Normally PP2A can remove a phosphate 
group from the threonine in ERK, thus causing a mixing of the layers 
but the mutation ensures that a phosphate of this kind is not present and 
PP2A cannot remove the phosphate from tyrosine. This leads to a cas¬ 
cade where the first layer allows two phosphorylation steps but the second 
only allows one. This cascade is represented in the following diagram: 


KKK 


KKK 



(90) 


KK KKP KKPP 



KKP'ase KKP'aie 


KKPP 



K 


KP 



Suppose that the phosphorylation and dephosphorylation in the first layer 
are distributive and sequential. In other words, only one phosphate group 
is added or removed in each reaction, the groups are added in a specified 
order and removed in the reverse order. Given this data it is possible to 
set up an MM-MA model as done in other cases above but this is not 
the model used in [25]. There the phosphorylation of MEK is modelled 
using mass action (MA) kinetics with both phosphates being added in 
one step and the concentration of Raf not included as a variable. On the 
other hand the action of MEK as an enzyme in phosphorylating ERK 
is modelled in detail. This gives a kind of hybrid MA/MM-MA model, 
which we call the PGS model. It is proved in that for the PGS model 
all solutions converge to a stationary solution at late times. 
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Now the MM-MA model for the in vitro system of [25] will be exam¬ 
ined, together with its MM reduction. In the Huang-Ferrell system discard 
the hrst four equations and take KKK to be a constant in the remaining 
equations. Then set [KP ■ KKPP], [KPP] and [KPP ■ KP'ase] to zero 
together with the reaction constant ag. Discard the last three equations. 
In the corresponding MM system this means taking KKK to be constant, 
setting [KPP] and 07 to zero and discarding the hrst and last equations. 
In this case the quantity KK is dehned to be [KKPP] + [K ■ KKPP], 
It is no longer possible to keep the coefficients with 7 < i < 10 

equal to the same positive constant 62 . This can be required for i 9 but 
K^^q must be replaced by zero. This has the effect that the expression 
b 2 {Ktot — [KPP]) is replaced in equations (iTOll . (ITTIl and d72ll by b 2 [K]. 
The system of equations obtained is 


_d ^ aiKKK[KK] 
dV l + bi{KKtot-Tn<) 

a2{KKtot-[KK]-Tn<) 

1 + 61 [kK,,, - [KK] - ’ 

d ™ aaliKKiKKtot - [KK] - TOi) 
dK l + bi{KKt^t-[KK]) 

aJ{K 

~b{KK -b (1 -b biiKKtot - [KK] - TO^)))! + 62 [if]) ’ 

d a5la<[K] 

dt^ ' l + b2[K] 
aeiKtot - [K]) 
l + b2{Ktot - [K])' 


(91) 


(92) 


(93) 


Next, in analogy with what has been done for other models above, sta¬ 
tionary solutions will be considered which satisfy the restriction [KK] = 
[KKPP]. These satisfy ai[KK]KKK = a 2 [KKP] and a 3 [KKP]KKK = 
a 4 [KK]. Hence aia^KKK^ = 0204 and KKK = J Substi¬ 


tuting this back in gives [KKP] 
^ /£i£ 4 ^ Hence 

V “2113 


aia 2 ^KKK[KK] = qi[KK] where 


A'A'tot = [KK] + [KKP] P KK = {2 + qi+b 2 [K]) [KK]. (94) 


This allows [KK] and [KKP] to be expressed in terms of KKtot, qi, fe 
and [K], These relations are equivalent to the hrst and second equations 
for stationary solutions. The remaining equation can be written as 

a,{K,,t-[K]){l+b2[K]) 

a5[K]{l + b2{Ktot-[K]))' ^ ’ 

One way of determinining a set of stationary solutions is as follows. First 
choose [K], Ktot > [K], & 2 , and the Oi. Then use the last equation to 
determine KK. Next use 


A'A'tot 


{2 + qT.+b 2 [K])KK 
l + b2[K] 


(96) 
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Then \KK], \KKP] and [KKPP] are determined in such a way that 
all the equations for stationary solutions are satisfied. To summurize, 
we get equilibria for parameters, conserved quantities and concentrations 
parametrized over [K], Ktot > [K], & 2 , and the ai by the following: 


KKK = 

/ a2a4 

V 

= 

a^{Kt.t-[K\){2 + q^+b2[K\) 

a5[K]{l+h2{Ktot-[K])) 

[KK] = 

a6(A'tot - [K]) 

as[K]{l + b 2 {Ki,t - [K])) 
a6{Ktot-[K]){l + b2[K]) 

lOi = 

a 5 [K]{l + b 2 {Ktot-[K]))’ 


(97) 

(98) 

(99) 

( 100 ) 


where qi = /£i£4 ^ 

In order to prove the existence of a Hopf bifurcation it is tempting to 
proceed as in the analysis of the MM system for the truncated MAP kinase 
cascade. To do so we need to find parameters such that the coefficients 
Aq, Ai and A 2 of the characteristic polynomial of the linearization are 
all positive. Then the parameters ai,i £ {1,2, 3,4} should be rescaled 
with a positive constant L so as to find a value Lq for which the Hurwitz 
quantity H 3 = A1A2 — Aq becomes zero. The linearization is very similar 
to that for the truncated MAP kinase cascade. Here again one entry is 
zero and all others have a sign. After a suitable change in the order of 
the variables only one sign differs. The resulting matrix of signs is 


- 0 
+ - 
+ - 


( 101 ) 


The main problem we have encountered in trying to implement this 
strategy is to find a point in parameter space where the positivity of 
the Ai holds. The coefficient Ao is the negative of the determinant of 
the linearization, while the coefficient A 2 is the determinant of the ana¬ 
logue of the submatrix M 2 defined in the case of the truncated MAPK 
cascade (using the modified order of the variables). Hence the Routh- 
Hurwitz method will give us the desired factor Lq only if the determinant 
of the linearization and that of M 2 have the same (negative) sign. Ex¬ 
periments with Maple indicate that these two determinants tend to have 
different signs for parameters satisfying the biologically motivated pos¬ 
itivity conditions. The signs are governed by the signs of polynomials. 
When attempting to attain the relevant combination of signs by fixing 
some of the parameters and varying others the signs of the determinants 
are governed by those of two polynomials. In all experiments we did these 
polynomials were different but shared a unique positive zero where their 
signs changed. Unfortunately, having the sign of one of the determinants 
negative requires being on one side of this zero while having the other 
negative requires being on the other side of it. After trying many param¬ 
eter configurations we are tempted to to conjecture that there is a deep 
reason why these determinants have systematically different signs. When 
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we apply the Routh-Hurwitz method in the case that Aq and A 2 have 
different signs it is still possible to find a unique positive value Lq of L 
with Hz = 0. However instead of leading to the purely imaginary eigen¬ 
values required for a Hopf bifurcation this leads to two real eigenvalues 
with equal magnitude and opposite sign. 


7 Conclusions and outlook 

The main result of this paper is a rigorous analytical proof of the existence 
of periodic solutions of the Huang-Ferrell system modelling the MAPK 
cascade. Their presence had been suggested by numerical and heuristic 
work in |26| . It was also proved that solutions of this type exist for 
cascades consisting of a layer with one phosphorylation followed by a layer 
with two phosphorylations but for the superficially similar case of a layer 
with two phosphorylations followed by a layer with one phosphorylation 
an attempt to obtain a similar proof ran into difficulties. The first of these 
two cases was considered in [2^ but to the authors knowledge the second 
was not previously investigated in the literature. The methods used in the 
proofs are bifurcation theory and geometric singular perturbation theory. 
It should be noted that the heuristic considerations in | 26| . which might 
in principle have been used as a basis for proofs, were in the end hardly 
used at all. It would be interesting to know whether this alternative route 
could also be effective in this problem. Relevant ideas, involving relaxation 
oscillations and the Conley index, are explained in [2], [12] and |13| . 

An important question is to what extent the mathematical results ob¬ 
tained here apply to the real biological system. First it should be noted 
that in nature there is not just one MAPK cascade but many. The ba¬ 
sic pattern is always the same but the details may be different. We now 
concentrate on the most famous example, the Raf-MEK-ERK cascade. 
This is essentially the example originally considered in HZ). In that case 
Raf is replaced by Mos but the architecture of the system is identical. In 
this model it is assumed that the phosphorylations and dephosphoryla¬ 
tions are distributive and sequential. Huang and Ferrell concentrate on 
this case but do mention that they also did simulations for the alterna¬ 
tive versions where one or both of the kinases act in a distributive way. 
In [26], |35| and the bulk of the present paper only the distributive and 
sequential case is considered. It has been found that phosphorylation of 
ERK by MEK is distributive but not sequential m, [3] while despho- 
sphorylation of ERK has been found to be distributive and sequential 
|34 |. On the other hand the phosphorylation of MEK by Raf has been 
found to be processive [I]. See also m where it is remarked that the 
distinction between the two mechanisms may not be absolute - processive 
phosphorylation may be thought of as a limiting case of distributive phos¬ 
phorylation where the second step takes place much faster than the first. 
This means that the original Huang-Ferrell model is not applicable to the 
Raf-MEK-ERK cascade. It is also known that there are cases in which 
processive phosphorylation suppresses complicated dynamical behaviour 
(in this case bistability) present when the phosphorylation is distributive 
m, 0. Thus it would be interesting to investigate whether there are os- 


22 


dilations in the system obtained by modifying the Huang-Ferrell model 
by making the phosphorylation in the second layer processive. 

It is also interesting to know what effect further interactions between 
proteins not included in the Huang-Ferrell model might have on the dy¬ 
namics. It was observed in m that binding of Raf to MEK can influence 
the dynamics of the MARK cascade, enhancing bistability. In real bio¬ 
logical systems the MARK cascade is also embedded in various external 
feedback loops. One well-known example is that ERK has a suppressive 
effect on Raf via the guanine nucleotide exchange factor son of sevenless 
(SOS). That the resulting negative feedback could lead to oscillations was 
observed theoretically in m- Sustained oscillations in the MARK cascade 
have been observed experimentally in [28]. They have a period of about 
15 minutes and have been observed to continue for over ten hours. A 
quantitative comparison with simulations indicates that these oscillations 
are not due to sequestration effects intrinsic to the cascade but to the 
feedback loop via SOS. Another type of feedback leading to oscillations 
which involves sequestration but is not intrinsic to the cascade is discussed 
in [^. In that case the binding of activated ERK to a substrate reduces 
its availability for dephosphorylation. These examples make it clear that 
there are numerous examples of biological interest which represent poten¬ 
tial applications of the methods developed in the present paper. 

Another key issue is that of the biological role of complicated dy¬ 
namical features such as bistability, sustained oscillations or chaos in the 
MARK cascade, with or without external feedback. Many different sig¬ 
nals pass through this cascade and it may be that non-trivial dynamics 
can be used to encode information, for instance by frequency modulation 
of oscillations. Here it could be useful to compare with other biological 
systems where this type of phenomenon is believed to be important, such 
as the NFkB pathway [55] or calcium signalling [7]. On the other hand it 
could be that complicated dynamical behaviour in the basic MARK cas¬ 
cade is an unwanted side effect and that the feedback loops in which it is 
embedded in biological systems serve to suppress it. Further mathemat¬ 
ical investigations of systems related to the MARK cascade could serve 
to understand the cascade itself better in its biological context and could 
also produce new insights into the architecture of biochemical systems. It 
should also be kept in mind that a better understanding of the dynamics 
of the cascade could be important for medical progress |31| . In chemother¬ 
apy of cancer Raf inhibitors are already in use while MEK inhibitors have 
been the subject of extensive clinical trials but have not yet been effective. 
A better theory of the system could help to understand where to look for 
appropriate drugs. 

One question which has been left open here is that of the stability 
of the periodic solutions whose existence has been proved. Is it possi¬ 
ble to develop methods to prove that in some of these models there are 
parameters for which the first Lyapunov coefficient is non-zero, which 
would solve the stability problem? Is it possible to extend the techniques 
used here to prove the existence of fold-Hopf or Hopf-Hopf bifurcations 
in the Huang-Ferrell model and to check the associated genericity con¬ 
ditions which would give information on chaotic behaviour? Evidently, 
modelling the MARK cascade gives rise to a large number of challenging 
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mathematical problems. 


Appendix: geometric singular perturba¬ 
tion 

theory (GSPT) 

In this appendix some results concerning GSPT needed in the paper will 
be collected. In m a theorem from [S] was applied but in this paper we 
need a parameter-dependent version of that result which does not obvi¬ 
ously follow from the theorem. In fact, starting from the basic transfor¬ 
mations carried out in and the idea of a slow manifold, the statements 
we need can be proved using standard results from the theory of centre 
manifolds. This will now be explained. 

The starting point is a system of equations of the form 

i = f{x,y,a,e), (102) 

<^y = g{x,y,a,t), (103) 

where x, y and a belong to open neighbourhoods of the origin in lR"i, 
and R*’ respectively and e belongs to an interval of the form [0, eo). 
The dot stands for d/dt. It will be assumed that the functions / and 
g are smooth and that they can be extended to smooth functions in a 
neighbourhood of e = 0. This system is now transformed as in m by 
defining a new time coordinate by r = t/e and treating e and a as new 
unknowns. The result is 


x'= ef{x,y,a,e), (104) 

y'= g{x,y,a,e), (105) 

a' = 0, (106) 

e' = 0, (107) 


where the prime denotes d/dr. Observe now that the solutions of the equa¬ 
tion g(x,y,a,Q) = 0 are stationary solutions of the system (I104D - (I107D . 
We assume that y(0, 0, 0,0) = 0. Suppose that there exists a smooth 
function ho such that g(x,y,a,0) = 0 is equivalent to j/ = ho(x,a). The 
centre subspace of the stationary point at the origin is of dimension at 
least ni -|- A: -I- 1. We now ensure that its dimension is no greater than that 
by assuming that all eigenvalues of the linearization of the system at the 
origin other than the zero eigenvalues arising from the manifold of sta¬ 
tionary solutions already mentioned and that coming from equation (fWl) 
have non-zero real parts. These will be called the transverse eigenvalues. 
They are the eigenvalues of Dyg(0, 0, 0, 0). For any positive integer I the 
centre manifold theorem ([19], Theorem 5.1) implies that there exists a 
centre manifold Me of the origin of class C* and of dimension ni + A: -|- 1. 
Another well-known result about centre manifolds states that any station¬ 
ary solutions sufficiently close to a given stationary solution must lie on its 
centre manifold. Thus the solutions of g{x, y, a, 0) = 0 all lie on the centre 
manifold of the origin. Since the dimension of the centre subspace is the 
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same for all of these points it follows that Me is also a centre manifold of 
these neighbouring points. Me is what is called the slow manifold. In a 
neighbourhood of the origin it can be written in the form y = h{x, a, e) for 
a C* function h with h{x, a, 0) = ho{x, a). Considering the restrictions of 
the dynamical system with the intersections of Me with the subspaces of 
constant a and e gives rise to a dynamical system depending in a regular 
way on the parameters a and e. Its explicit form (when written in terms 
of the time coordinate t) is 

X = f{x,h{x,a,€),a,e). (108) 

In this way the singular limit e —>■ 0 in the original system has been 
reduced to a regular limit. For e = 0 it reduces to 

X = f{x,ho{x,a),a,0). (109) 
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